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This paper discusses the formulation and development of a trajectory reconstruction 
tool for the NASA X— 43A/Hyper— X high speed research vehicle, and its implementation 
for the reconstruction and analysis of flight test data. Extended Kalman filtering techniques 
are employed to reconstruct the trajectory of the vehicle, based upon numerical integration 
of inertial measurement data along with redundant measurements of the vehicle state. The 
equations of motion are formulated in order to include the effects of several systematic 
error sources, whose values may also be estimated by the filtering routines. Additionally, 
smoothing algorithms have been implemented in which the final value of the state (or 
an augmented state that includes other systematic error parameters to be estimated) and 
covariance are propagated back to the initial time to generate the best-estimated trajectory, 
based upon all available data. The methods are applied to the problem of reconstructing 
the trajectory of the Hyper-X vehicle from flight data. 


I. Introduction 

T he NASA Hyper-X Program is a joint Langley Research Center and Dryden Flight Research Center 
effort to develop a hypersonic airbreathing research vehicle, powered by an airframe-integrated supersonic 
combustion ramjet (scramjet) propulsion system. Major objectives of the program are to acquire data on the 
operation of the scramjet propulsion system and to demonstrate controlled hypersonic flight technologies. 1,2 
A secondary goal is to evaluate the experimental 3-6 and computational' -10 techniques used in the design and 
analysis of the Hyper-X vehicle in order to facilitate future hypersonic and launch vehicle programs. 1,11,12 

The flight experiment mission profile is shown in Fig. 1, and takes place in the Western Test Range over 
the Pacific Ocean. The experiment begins with a boost phase, where the Hyper-X Launch Vehicle (HXLV) 
stack (shown in Fig. 2), consisting of the X-43A/Hyper-X Research Vehicle (HXRV), shown in Fig. 3, mated 
to a modified Pegasus booster, is dropped from the NASA B-52 aircraft at an altitude of approximately 
40,000 ft and ascends to an altitude of approximately 100,000 ft. At this point, the HXRV separates from 
the stack, stabilizes, and the scramjet propulsion experiment begins, lasting several seconds. Following the 
engine test, the vehicle follows a descent trajectory that includes various maneuvers that are optimized for 
the post- flight extraction of stability and control derivatives and other aerodynamic properties (e.g., Ref. 13). 
It is important to point out that approximately one minute prior to the drop time, the on-board navigation 
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Figure 1. Hyper— X Flight Trajectory 



Figure 2. Hyper— X Launch Vehicle Stack Configuration 


system switches from a satellite-aided system to a free inertial system. Therefore, drift is expected in the 
on-board navigation solution after this time. 

A critical component of the post- flight data reduction and analysis tasks to determine propulsion sys- 
tem performance and the extraction of various aerodynamic coefficients is the reconstruction of the vehicle 
trajectory. The trajectory reconstruction process uses in- flight telemetry data combined with a variety of 
redundant measurements to generate maximum-likelihood estimates of the vehicle position, velocity, and 
orientation time histories, along with estimates of the uncertainty associated with them. The Hyper-X vehi- 
cle has many instruments that gather data to be used for this process, such as the Inertial Measurement Unit 
(IMU), Global Positioning System (GPS), and Flush Air Data Systems (FADS). Numerous other useful ob- 
servations and measurements are available from ground-based radar tracking stations and from instruments 
carried on board the B-52 carrier vehicle. 

A tool known as the Statistical Trajectory Estimation Program (STEP) 14-16 has been developed for 
post-flight trajectory reconstruction and parameter identification in support of maneuvering re-entry vehi- 
cle flight projects in the 1960’s. STEP uses a recursive, minimum variance filtering algorithm to process 
accelerometer and gyroscope measurements recorded by the IMU during the flight, as well as position and 
velocity measurements from outside sources, such as ground based radar tracking stations, to determine 
a maximum-likelihood trajectory and user defined systematic error sources along with an estimate of the 
accuracy of these parameters. STEP allows the user to specify which systematic error sources should be 
estimated, such as an accelerometer bias, and is general with the type of outside tracking measurements 
available. STEP and its variants have been applied to a variety of post-flight planetary entry trajectory re- 
construction and analysis problems, particularly the SV-5D/X-23A maneuverable entry vehicle, 17 the Viking 
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Figure 3. Hyper— X Research Vehicle Configuration 


landers, 18-20 Pioneer, 21 and the Space Shuttle, 22 ’ 23 among others. 

STEP is not without drawbacks, however. The STEP user interface is poor, and the source code is 
archaic and not commented. The latter issue makes the STEP algorithm difficult to check or modify to 
suit algorithmic developments or new types of data sources. With these limitations in mind, an effort 
was made to assess how STEP may be modernized both in terms of its underlying filter algorithm and 
its implementation. 24 This paper discusses the application of several nonlinear filtering algorithms to the 
post- flight trajectory reconstruction problem, wherein the inertial accelerations and angular rates measured 
by the IMU are used to integrate the nonlinear equations of motion, while GPS observations of position and 
velocity are used to provide corrections to the trajectory estimate. These corrections are incorporated into 
the state estimate at the measurement points through use of an Extended Kalman Filter (EKF) 25-28 for a 
forward and backward pass through the data. The results of the two filtering passes are combined using the 
Fraser-Potter smoothing method 26,29 to determine the final trajectory estimate as a function of time. 

Similar methods have been recently applied to the entry trajectory estimation problem in Refs. 30-32 
using a variety of Kalman filtering and smoothing techniques. Unfortunately, these methods do not model 
the vehicle rotational dynamics, and furthermore, Refs. 31 and 32 are not formulated for cases where the 
vehicle measures its acceleration and angular rates. Therefore, the methods of Refs. 30-32 are not directly 
applicable to the problem investigated in this paper. Other proposed techniques for post-flight trajectory 
reconstruction simply use either a numerical integration of the raw IMU acceleration and angular rate data, 
or the on-board navigation solution, to produce estimates of position, velocity, and orientation, 33, 34 but do 
not make use of all available observation data through the application of a filtering and smoothing process 
to correct for various drifts, as is discussed in the next section. 

II. Filtering Equations 

The current section shall be concerned with the problem of estimating the state of a nonlinear system of 
ordinary differential equations 

x(i) = f(x(t),v(t),t) (1) 

where x is the matrix of the system states, f(x,v,f) is a nonlinear function that describes the evolution of 
the system state, and v is a random process with covariance matrix Q given by Q = i?(vv T ) and mean 
v = -E(v), where E is the expectation of the random variable. The state estimation procedure requires 
knowledge of independent observations of some quantity which may be related to the value of the state by 
the nonlinear algebraic equations 

yW = g(x(f),w(t),f) (2) 

given at discrete time intervals, where y is the value of the observation, g is the nonlinear relationship 
between the states and the observations, and w is a random measurement noise variable with covariance 
matrix R and mean w. 
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A. The Extended Kalman Filter 


The extended Kalman filter (EKF) is a somewhat ad-hoc modification of the linear Kalman filter in which the 
linear state transition is replaced by a numerical integration of the full nonlinear dynamic model in between 
measurement updates, before applying the Kalman update at each measurement time. In addition, the 
nonlinear relationship between the state and measurement data is used to calculate the residuals rather than 
the linear approximate relationship. The EKF uses the current estimated state as the reference conditions 
for computing the system dynamics matrices required to calculate the Kalman gain matrix at each time 
step. The algorithm is given by 25-28 


Xfc = x fe + K fc [y fe - g (x fe ,w fe ,ffc)] (3) 

where x^- is the predicted value of the state at time increment k, y k is the value of the measurement data, 
and g (xfc, Wfc,ffc) is the predicted value of the measurement data at time t k , following from the propagation 
of the state estimate from time t k -\ to time t k . The gain and covariance equations are given by 

K k = PfcCjT (CfePfcC^ + D fc R fc D^) _1 (4) 

P k = (I — KfcCfc) Pfc (5) 

where is the solution of the matrix differential equation 

P(i) = AP + P t A t + BQB t (6) 

P(tfc-r) = Pfc-i (7) 


at time t k , calculated through numerical integration. The matrices A, B, C and D used in the Kalman gain 
calculations are found by linearizing the system dynamics about the current state estimate, x/ ;; , so that 


A(t) 
B (t) 
C(t) 
D(t) 


df_ 

9x 

df_ 

dv 

dg 

9x 

dg 

dw 


x=Xfc ,v=0 


x=Xfc , v=0,i 


X-Xl; ,W = 0 


x— Xfc ./w—O.t 


(8) 

(9) 

(10) 

( 11 ) 


The EKF formulation is superior to that implemented in the original STEP code in Refs. 14-16. The 
STEP implementation made use of linearized equations of motion for state propagation between measurement 
points, as well as linearized representations of the measurement equations for computing the residuals. This 
linearization was performed about the current estimate, as is used by the EKF for gain calculations, or 
about a pre-defined nominal trajectory; the latter technique is better known today as a neighboring-optimal 
Kalman filter 2 ' or as a linearized Kalman filter. 28 

Note that in the current formulation of the filtering equations, the state x could include both dynamic 
states as well as constant model parameters to be estimated. Likewise, the input noise v could include 
random inputs to reflect the model uncertainty (process noise), or the uncertainty associated with model 
parameters that are not included in the state x. The uncertain parameters from the latter case are commonly 
known as “consider parameters.” The two noise inputs to the system dynamics both seek to degrade the 
state covariance matrix in order to reflect uncertainty in the model, but are different in the sense that the 
process noise accounts for effects that are unknown whereas the consider parameters account for effects that 
are known but uncertain. A Kalman filter designed to include parametric uncertainty is typically known 
as a Kalman-Schmidt filter. 25 It is worth noting that the STEP filter was designed to have the option of 
using a variety of consider parameters, but did not include a means of incorporating process noise into the 
estimation procedure. 

The linearized Kalman-Schmidt filter algorithm implemented in STEP is computationally efficient and 
well understood. However, the algorithm has proven itself to be sensitive to the quality of the nominal 
trajectory or the estimated trajectory, depending on the mode of operation. If the reference trajectory is not 
an accurate representation of the true flight path, then the assumption of local linearity is violated and the 


4 of 21 


American Institute of Aeronautics and Astronautics 



filter often diverges, hence the algorithm must be modified to take nonlinearity into consideration, which is 
the advantage of the EKF method. 

The EKF algorithm implemented in this paper may be used for both a forward pass and a backward 
pass through the measurement data. The final estimate from the forward pass through the data is used as 
the initial condition for the backward pass. This process can also be iterated globally using the change in 
the initial condition to measure the convergence of the estimates. The final trajectory is then computed as 
a linear combination of the forward and backward passes using the Fraser-Potter method, 26, 29 which is the 
subject of the next section. 


B. The Fraser Potter Smoother 


The optimal smoothed solution is the blend of the forward and backward estimates of the state time history as 
output from the EKF algorithm. The Fraser-Potter method 26, 29 seeks to determine the smoothed solution 
at any instant in time as a linear combination of the forward and and backward filtered solution. The 
smoothed solution at a time k is given by 


Pfe = 



(12) 


Xfc 


Pk 




(13) 


where P is the smoothed covariance estimate, and x is the smoothed state estimate, P/ and P^ are the 
covariance estimates from the forward and backward filter, respectively, and x.f and xj, are the state estimates 
from the forward and backward filter. By inspecting Eqs. (12-13) it is apparent that the smoother combines 
the forward estimate at time tk, which makes use of all measurement data from the beginning of the data 
set to time tk , with the backward estimate at time tk, which makes use of all data from time tk to the end 
of the interval. In other words, the smoother makes use of all available data at each point of interest. 

It should be noted that the original STEP code did not make use of an statistical smoothing method as 
is discussed in this section. STEP conducted smoothing by simply integrating the final state estimate back 
to the initial time via Eq. (1). The flaw with this approach is that it does not account for model uncertainty, 
and is also subject to numerical error buildup over large fit spans. These problems are eliminated by using 
the forward backward filter implementation, combined together in the smoother. 


III. Dynamic Model 

The dynamic model discussed in this section is similar to that implemented in the original STEP. The 
equations of motion are based upon Newton’s laws of motion, and are developed in a local, non-inertial 
coordinate system. Similar to the formulation used in STEP, the inertial acceleration and angular rate 
measurements recorded during the flight are used in place of detailed aerodynamic and propulsion force and 
moment models. An approximate representation of the complete gravitational potential is used to develop 
an expression for the gravitational force acting upon the vehicle. These assumptions lead to a system of 
nonlinear ordinary differential equations that describe the motion of the vehicle. The summary to follow 
begins with a definition of the required coordinate systems. 

A. Coordinate Systems 

This section provides a brief review of the coordinate systems used to develop the nonlinear equations of 
motion used in the trajectory reconstruction process. Firstly, the Earth Centered Inertial ( ECI ) frame is 
the frame fixed at the center of the Earth, or any other primary body of interest. The fundamental plane 
of this coordinate system is the equatorial plane of the Earth, with the x-axis, Bxi, aligned with the vernal 
equinox, the 2-axis, bzi, in the direction of the planetary axis of rotation, and y- axis, Byi, chosen in such a 
way as to complete a right-handed coordinate system. This frame is obviously not truly inertial, although 
for short periods of time the inertial assumption is a valid and common approximation. 

The Earth Centered Earth Fixed ( ECEF ) frame is the coordinate system whose origin is located at the 
center of the Earth, similar to the ECI frame, but the coordinate axes are fixed with respect to the planet 
surface. The x-axis, bxf, and the 2 -axis, b Z f, are aligned with the prime meridian and the axis of planetary 
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rotation, respectively. The ECEF and ECI frames may be related by the simple transformation 


( > 


f } 

&XF 


exi j 

1 - 

I l YF 

> = GECI2ECEF < 

- 1 

e YI 

l ® ZF , 


, j 


GeCI2ECEF 


cos A sin A 0 
— sin A cos A 0 
0 0 1 


(14) 


(15) 


where A is the sidereal angle of the prime meridian. 

The geocentric ( GC ) frame is one whose origin is located at the center of mass of the vehicle, and whose 
2 -axis, ezc , is aligned in the direction of the origin of the ECI (and ECEF) frame. The cc-axis of the this 
frame, &xCi is perpendicular to the 2 -axis, and points in the general direction of north. The GC frame may 
be related to the ECEF frame by the transformation 
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1 - 
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> = GeCEF2GC ' 

-> 1 
l YF 



, ® ZF J 



— sin (j> 

0 

COS (j) 


GeCEF2GC = 

0 

1 

0 



— COS (j) 

0 

— sin (f> 



cos 9 
— sin 9 
0 


sin 9 
cos 9 
0 


0 

0 

1 


(16) 


(17) 


where 9 is the longitude of the vehicle and 4> is its declination. The geocentric frame is sometimes called the 
topocentric or local-vertical frame. The ECI, ECEF, and GC frames are shown in Fig. 4. 

A frame similar to the GC frame is the geodetic (CD) frame, sometimes known as the topodetic, local- 
horizontal, or North-Earth-Down (NED) frame. This frame is one in which the origin is located at the center 
of mass of the vehicle, and the rr-axis, points in the direction tangent to the surface of the Earth, in 
the direction of North. The 2 -axis, §zDi is aligned perpendicular to the araxis and oriented downward, 
while the y- axis, §yd is shared with the GC frame. The CD frame may be related to the GC frame by the 
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Figure 5. Declination, Geocentric Latitude, Geodetic Latitude, Geodetic Frame, and Geocentric Frame 


transformation 


f C,v D 


exc 1 

1 - 

1 Gyd 

> = GgC2GD < 

- 1 

e YC 

{ e ZD j 


, ^zc J 


GgC2GD 


cos e 0 sin e 
0 10 
— sin e 0 cos e 


(18) 


(19) 


where e = (f> g d — <p, with <j> g( i being the geodetic latitude of the vehicle’s sub-latitude point, which is the 
projection of the vehicle position to the surface of the Earth along the e ZD axis. The geodetic latitude is 
then the angle between the e Z D axis and the equatorial plane of the Earth. The GC and GD frames are 
shown in Fig. 5. 

An alternate measure of the location of the sub-latitude point is the geocentric latitude, <f> gc , which is the 
angle from the equatorial plane of the Earth to the line connecting the sub-latitude point and the center of 
the Earth. Fig. 5 shows the differences between the declination, geocentric latitude, and geodetic latitude 
for a planet with exaggerated eccentricity. Note that the three angles are equal for a spherical planet. Fig. 5 
also shows the physical meaning of the geodetic altitude, h gr j, and the magnitude of the radius vector, r - 
two variables that will be used presently. 


Table 1. Summary of WGS— 84 Earth Model Constants 


Constant 

Value 

-Re 

2.092565 -10 7 ft 

e e 

8.1819191 -10" 2 


1.4076444 -10 16 ft 3 /s 2 

J 2 

1.0826230 -10" 3 

0 

4.1780741 -10" 3 °/s 


The relationship between the geocentric and geodetic latitude is quite simple, and takes the form 

tan 4> gc — (l ^ 0 ) tan (f) g d (20) 
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where e 0 is the Earth’s eccentricity. The transformation between declination and geodetic latitude, however, 
is more complicated. Vallado 35 gives the formulas 


tan 4> gd 


hgd 


tan (j> 


sin <j>gd 


r cos 4> ijl - e\ sin 2 (j> gd 
r cos 4> f?® 

C0S ^ sd ^/l — e® sin 2 ct>g d 


(21) 

(22) 


where i?® is the mean equatorial radius of the planet. Sofair 36,37 determined an exact solution of Eqs. (21- 
22), (f> g d (r, (j> ) and h gd (r, 0). The solution is explicit, and is free of singularities at the poles and the equator, 
although the solution space is confined to a region outside a spheroid fixed at the center of the Earth with 
semimajor axis of approximately 1.42 x 10 5 ft. This restriction is of no concern for the present problem as 
the region is well within the surface of the Earth. 

The transformation to geocentric coordinates is much less complicated, and may be determined directly 
from Eqs. (21) and (22). Note that there are other definitions of the geocentric coordinates available in the 
literature. For an example, see Ref. 38. 

The values of i?®, e®, and other constants that describe the shape, mass, and rotation of the Earth are 
summarized in Table 1. These values are the World Geodetic System of 1984 (WGS-84). 39 

Lastly, the body ( B ) frame is fixed to the center of mass of the vehicle with the §xb axis directed through 
the nose of the vehicle, the §zb axis oriented perpendicular to the §xb axis, in the plane of symmetry, with 
the positive direction taken to be toward the bottom of the vehicle, and the eys axis completing the right- 
handed system so that it points out the right side of the vehicle. The B frame is shown in Fig. 3. The 
orientation of the B frame is typically expressed by a 3-2-1 Euler angle sequence relating the GD frame to 
the B frame. The transformation is given by 


f e X B 


&xd 1 
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evs 

> = G GD2B * 

3 

l ® ZB / 


e zd j 


’ 1 

0 

0 


COS0 

0 

— sin 0 


cosT 

sin T 

0 ’ 

0 

cos<3> 

sin <J> 


0 

1 

0 


— sin T 

cos T 

0 

0 

— sin $ 

cos $ 


sin 0 

0 

COS0 


0 

0 

1 


(23) 


(24) 


where T, 0, and $ are the azimuth, elevation, and bank angles, respectively, also commonly referred to as 
the yaw, pitch, and roll angles. 

Although the Euler angle sequence offers an ease of visualization, other parameterizations of the vehicle 
orientation are preferable for implementation in the equations of motion and filtering algorithms. 40 The most 
advantageous of these parameterizations are the Euler parameters, the use of which leads to a representation 
of the attitude transformation matrix of the form 


GqD2B — 


P 2 1 P 2 _ 2 _ 2 

e 0 w c 1 e 2 e 3 
2 (eie2 — eoe 3 ) 

2 (eie 3 + eoe2) 


2 (eie 2 + e 0 e 3 ) 

•0 G t e 2 e 3 
2 (e2e 3 — eoei) 


2 (eie3 — eoe2) 

2 (eoei +e 2 e 3 ) 

^2 „2 _i_ „2 

e 0 C 1 c 2 ' c 3 


(25) 


The Euler parameters constitute a four parameter set that describe the orientation of the B frame with 
respect to the GD frame. Due to the fact that four parameters are used to describe a rotation with three 
degrees of freedom, the constraint equation e 3 + e\ + e 2 + e§ = 1 must also be satisfied. 


B. Equations of Motion 


Newton’s laws of motion, expressed for a body moving relative to a planet- fixed reference frame, lead to the 
translational equations of motion of the form 24 


f 


9 


—w 

u 


r 

v 

rcoscj) 


- n 


(26) 

(27) 

(28) 
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V 

w 


1 / 2 3 f-L J 2 , , i . 

axe H — (mw — v tan 0) j- sin (20) 

arc H — (uv tan 0 + vw) 


1 / 2 , 2\ , P 

aze (« + f ) + 


3/iJ 2 

2r 4 


(2 — 3 cos 2 0) 


(29) 

(30) 

(31) 


where it, u, and u; are the x, y , and z-axis components of the velocity in the GC frame, respectively, and 
axe, o-yc, a-zc are the geocentric representation of the accelerations acting at the body center of mass, not 
including the gravitational terms. 

These accelerations may be given in terms of the body-axis accelerations at the IMU location by the 
transformations 


&xb 

ayB 

... 

ax t 
ay T 

► - (r 2 + i 1 ) < 

xm 

Dm 

> -2T < 

xm 

Vm 


xm 1 
Vm j 

azB J 


k «ZT y 


, Z M , 




, z m J 


(32) 


where axB, ays, and uzb are the accelerations acting at the vehicle center of mass, represented in the B 
frame, axr, ayr, and azr are the accelerations acting at the IMU location, expressed in a frame parallel 
to the B frame, Xm, l I m, and zm are the components of the position of the IMU measured from the vehicle 
center of mass, expressed in the B frame, and 


r = 


0 -R Q 
R 0 -P 
-Q P 0 


(33) 


where P, Q, and R and the components of the angular velocity of the B frame with respect to the ECI 
frame, expressed in the B frame. 

The accelerations may then be transformed into the GC frame to be consistent with the equations of 
motion by means of 

! axc I j a XB 

a YC ( = £*B2GC \ ay B 

aze ) [ «zb 

where G B 2 GC = (Ggd 2 bGgc 2 Gd) T ■ 

The rotational kinematics of a vehicle relative to the GD frame may be written as 




(35) 


In summary, the equations of motion for a vehicle relative to a rotating central body are given by the 
combination of Eqs. (26-31) and (35), together with the transformation given in Eq. (32). The states may 

be written as x = | r 0 9 u v w eg e\ e 2 j ■ The strapdown formulation of the inertial 
measurements in the equations of motion may be modified without difficulty to a platform implementation 
as in Ref. 41 for applications that make use of that particular type of inertial measurement source. 

The dynamic model discussed in this section is similar to that implemented in the original STEP code 
with several subtle differences. The first and most obvious difference is that the current implementation uses 
the radius of the vehicle position relative to the planet center of mass whereas the STEP equations of motion 
used an altitude variable. The current formulation is superior in the sense that additional transformations 
between altitude and radius as required to compute the gravitational acceleration is not required. Another 
difference between the current dynamic model and the STEP dynamic model is that the latter did not 
distinguish between the GD frame and the GC frame as does the former. Although the difference between 
these two frames is a small, it is not negligible, and may lead to large errors in the north and down velocity 
components, as was demonstrated in Ref. 42. 

It should also be noted that the equations of motion developed in this section make use the Euler 
parameters to describe the orientation of the P-frame relative to the GD-frame, as was used in STEP as 
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well. Recent research has shown that some problems exist when using this particular parameterization 
within estimation algorithms, which lead to the possibility of ill-conditioned covariance matrices due to the 
constraint on the norm of the Euler parameters. 43-45 Pittelkau, 46 however, showed that filters with properly 
tuned process noise models do not suffer from these problems, and that enforcing the Euler parameter 
constraint within the filter update equation is not an ad-hoc modification of the filter, but in fact equivalent 
to a measurement update equation. Therefore, the use of the Euler parameter set within the estimation 
procedure remains justified. 

The next section discusses the development of the equations relating the states to the measurement data 
for use in the filtering process. 

IV. Measurement Equations 

This section discusses the formulation of the measurement equations used within the trajectory recon- 
struction tool, including the modelling of various error sources. The first such measurement equation is that 
of the GPS, or some other source such as a radar altimeter, that provides redundant position information 
that may be incorporated into the trajectory estimate during the filtering process. Secondly, a GPS unit 
or other sources may provide a velocity measurement, which may also be incorporated into the trajectory 
estimation procedure. Another source of interest is the IMU data, which records the acceleration and angular 
rates of the vehicle in question. This data source in not used to provide measurement updates via a filtering 
algorithm, but are inputs into the equations of motion that are integrated to provide the state predictions at 
the time of the measurement update. However, a vehicle that contains two separate IMUs could use one to 
process the equations of motion, and the other to provide redundant data for use in filtering. Also included 
in this section are equations for the measurement of the vehicle orientation, and for ground based radar 
tracking stations. 

A. Position Measurement Equations 

The filtering algorithms require observations of the vehicle state in order to generate a minimum variance 
trajectory estimate. One such source of the observation data is a satellite navigation Global Positioning 
System (GPS) receiver carried within the vehicle. The GPS-based observations provide position estimates 
of the receiver location. 

The position observation data, provided in the geodetic coordinates h g d M , 4>gdMi and ^gd M i can be written 
as a function of the vehicle state by the transformations 

f x* | ( xg + Sxg 1 

N y* f = G B2 ECEF < Ug + Syg f 

[ v J [ zg + Szg ) 

x f g = r cos (j) cos 8 + x* 
y f g = r cos (j> sin 6 + y* 

Zf g = r sin <j>+z* 


r G 

= V x2 F G + y 2 F G + z k 

(40) 

sin tj>Q 

- (%) 

(41) 

cos 0q 

Xf g 

r G cos (j) G 

(42) 

sin 0q 

Vf g 

re cos <f>G 

(43) 

hgd M 

= h gd {r G ,(t>G ) 

(44) 

4 gd M 

= <t>gd(r G ,<l>G ) 

(45) 


where xq, yc, and zq are the components of the position of the GPS unit relative to the vehicle center 
of mass, expressed in the B frame (with bias error sources 6xg, Syg, and Szg), while x*, y *, and z* are 
the same position expressed in the ECEF frame. The variables xf g , Vf g , and zf g represent the ECEF 
coordinates of the GPS receiver, Fq, <j>G, and Oq are the geocentric coordinates of the GPS receiver, and 
h g d M and <j) g d, M a re the position measurements, calculated using Sofair’s method. 
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B. Velocity Measurement Equations 


Similarly, the velocity of the GPS unit or other observation is provided in the planet-relative, geodetic 
velocities u g d M , Vgd M , and w g d M , and may be found as a function of the inertial velocity of the vehicle center 
of mass, in the GC frame, by first transforming the inertial velocity from the GC frame to the GD frame, and 
then subtracting out the velocity due to the planetary rotation in order to create planet-relative velocities. 
Lastly, the velocity due to the rotation of the body and the offset of the GPS unit from the vehicle center of 
mass must also be accounted for. These transformations take the form 

I u gd,M 1 f u 1 f xq + Sxg 1 

\ v gd M f = G GC 2 GD < V — rtt COS (j) /+G_B 2 GDr< y G + S YG > (46) 

{ WgdM j { w J { zg + Szg J 


C. Acceleration and Angular Rate Measurement Equations 


The equations of motion for a vehicle, developed in Sec. Ill B, assume that the acceleration and angular 
rates of the vehicle center of mass are known functions of time. The inertial flight data contains this 
information, however, several transformations are required in order to bring this data into the proper form, 
and furthermore, the data recorded during the flight is subject to systematic errors and measurement noise. 
Systematic error sources may be modelled, and if values are known, their effects may be compensated for 
in the trajectory reconstruction procedure. Otherwise, the values of the systematic errors may be estimated 
using the filtering techniques, with an augmented state vector consisting of the usual position, velocity, and 
attitude states, along with the systematic error terms that are assumed to be constants over the time interval 
in question. 

Assuming small systematic error sources allows for a linear representation of the IMU measurement 
equations of the form 16, 47,48 


j axM 


| &YM 

> = 

{ a Z M J 



1 +£x 
Vyx 
fJ-zx 


1-iXY 

l-l XZ 


f axr j 


f VX } 

1 + 6- 

l-l YZ 

S 

1 1 

1 QYT 1 

[ + 1 

VY 

HZY 

i + iz 


[ &ZT J 


{ V Z ) 


(47) 


where axr, clyt , and azr are the true accelerations at the IMU location, expressed in the B frame, and 
a.Y m, <1ym , and azxi are the measured accelerations at the IMU location, expressed in a frame parallel to 
the B frame. The parameters 6f> 6^> and £z account for the scale factor error in the accelerometer data, 
Hxy, l^xz, Hyx, I^yz, [dzx, and yzY account for misalignment errors. The parameters ijx , t/y, and rjz 
account for bias. The components of the position of the IMU with respect to the vehicle center of mass (xm, 
UMi and Zm) are also assumed to be perturbed by the biases Sxm, Sym, and Szm, respectively. 

Similarly, the inertial angular velocity components are subject to systematic errors of the form 16,4 ' -48 


Pm j 


l + £p 

HPQ 

Qm | 

f = 

HQP 

l + 

Rm \ 


M RP 

fJ-RQ 


Upr 
HQR 
1 + 6? 



(48) 


where Pm , Qmi and Rm are the measured inertial angular velocity components, expressed in the B frame, 
and P, Q, and R are the true inertial angular velocity components of the vehicle, expressed in the B frame. 
£p, , and 6? are the scale factor errors, ypQ, ypp , /i qp , hqr, /i rp , hrq are the misalignment errors, and 

VPi V Qi and rjp are the bias terms. Note that the inertial angular velocity is a property of the body and not 
a point, and so these terms do not require a moment-arm transformation as do the accelerations. 

Note that the preceding transformations give the measured acceleration and angular rate data as a 
function of the true acceleration and angular rates, thus, in order to provide the proper inputs into the 
equations of motion, the IMU measurement equations must be inverted to yield the true acceleration as a 
function of the measured acceleration and error sources, either known or estimated. 

The acceleration measurement error equations developed in this section are used principally to correct 
the forcing function input to the equations of motion. However, a vehicle with two or more IMUs could 
potentially use the additional data as a source of redundancy, which may be incorporated into the state 
estimation in a manner similar to the position and velocity measurements. 
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D. Radar Measurement Equations 


The last measurement source of interest to the current document is that associated with ground-based radar 
tracking stations that measure the range, azimuth, and elevation of the vehicle with respect to a GD-frame 
local to the tracking station. Given the site coordinates h gds , 4>gds, Os, or equivalently, rs, 4>s and Os, the 
relative position between the vehicle and the tracking station, expressed in a GD-frame local to the station, 
may be expressed as 

I x re i I I r cos </> cos 0 — rs cos 4>s cos Os 1 

< Urei / = < r cos (j> sin 0 — rs cos (j>s sin Os / (49) 


Zrel 


r sin <t> — rs sin <j>s 


where 



- sin (j> gds 0 cos (j) gds 


cos Os sin 6s 0 

G s = 

0 10 


— sin Os cos Os 0 


-cos 4> gds 0 — sin 4> g d s 


0 0 1 


The range, azimuth, and elevation angles may then be written as 


(50) 


V 

sinx 

sin k 

cos K 


\J X rel + Orel + Z rel 

Zrel 

V 

Vrel 

V X r el + Orel 

%rel 

V X rel + y‘rel 


(51) 

(52) 

(53) 

(54) 


where v is the range measurement, \ is the elevation angle, and k is the azimuth angle, respectively. 


V. Application to the Analysis of Flight Data 

This section summarizes the results of applying the methods developed in the previous sections to the 
problem of estimating the trajectory of the Hyper-X flight test vehicle given the IMU data and several 
redundant data sources, such as GPS and radar measurements. The data preprocessing methods are discussed 
first. 


A. Flight Data Preprocessing 

This section summarizes the data preprocessing that were required in order to provide useful measurements 
as needed for the trajectory reconstruction process. The raw flight data were passed through a preprocessing 
tool, known as PREPR, 49 which consists of four main utilities. The first utility removes duplicate points in 
the data series which result from differences between the sensor data refresh rate and the telemetry sampling 
rate. The second utility allows users to interactively remove rogue points from the data series. Users select 
individual or groups of points to remove from a graph of the data series. The third utility samples the 
data series at a uniform time step, using a user-specified option of linear, cubic spline, or Hermite spline 
interpolation between available data points. The fourth utility gives the user the option to filter the data 
series using Fourier series and optimal Weiner filtering method 50 available within the Systems Identification 
Program for Aircraft (SIDPAC). 51 A separate fifth utility is included which can be used to differentiate a 
data series once it has been passed through the previous four routines. The user is given the option to use 
backward difference or center divided difference methods. Alternatively, users can select to curve fit the data 
using either cubic or Hermite spline polynomials, and to individually differentiate each piecewise segment. 

Data recorded by the Hyper-X control room were used for the drop phase up to the end of the engine 
test experiment. Some communications issues arose during the descent portion of the flight in which data 
were not transmitted to the control room, but were recorded by a P3 chase plane. The data from the control 
room and from the P3 chase plane were spliced together to give a continuous data stream for post-flight 
analysis. This spliced data set was used for the portion of the flight beginning just after the engine test to 
the end of the flight. 
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Figure 6. Body Y-Axis Inertial Angular Velocity 


1. Inertial Data 

Measurements of acceleration and angular rates were provided by QA2000 accelerometers and GG1320 Digital 
Laser Gyroscopes (DLG) contained within the Honeywell H-764G system, 52 a unit with a background of 
numerous civilian and military applications. These measurement data were obtained from the on-board 
telemetry stream at a rate of 100 Hz and were first preprocessed by removing points which contained 
obviously bad or rogue data. Next, the data were edited to remove points which could be considered as 
repeated, stale, or non-updated frames. Gaps in the data due to telemetry dropouts or due to edited data 
were filled in by a piecewise cubic Hermite interpolating function, and were then filtered to remove the effects 
of high-frequency structural noise, measurement noise, and to smooth out the data from quantization effects. 

Table 2. Summary of Cut-Off Frequencies for IMU Data Preprocessing 


IMU Data 

HXLV Stack Cut-Off 

HXRV Cut-Off 

a-XM 

10 Hz 

10 Hz 

O-YM 

8 Hz 

10 Hz 

a-ZM 

8 Hz 

10 Hz 

Pm 

8 Hz 

10 Hz 

Qm 

5 Hz 

10 Hz 

Pm 

8 Hz 

10 Hz 


Table 2 summarizes the cut-off frequencies used to filter the IMU data for the periods when the HXLV 
stack is complete and for the periods after the HXRV has separated from the HXLV. The body Y-Axis inertial 
angular velocity is shown in Fig. 6 as an example of the data before and after the filtering process, during 
the boost phase. In this example, evidence of a known 8 Hz structural mode in the HXLV configuration can 
clearly be seen, which influenced the decision to use a lower cut-off frequency for this channel. After the 
separation of the HXRV from the HXLV, the IMU data had almost no apparent structural mode content, 
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and therefore a higher cut-off frequency was used for preprocessing. 

Numerical derivatives of the inertial angular rates, as required for evaluation of Eq. (32), were calculated 
using a second-order centered divided difference scheme, except for the beginning and end points, which 
used first-order forward and backward divided differences, respectively. 

Preliminary values of the process noise standard deviations used in connection with the inertial measure- 
ment data in the EKF formulation were taken from Refs. 53 and 54, which were tuned somewhat to suit the 
problem at hand. 

2. GPS Data 

GPS measurements were obtained as outputs from the H-764G GPS Embedded Module (GEM), which is 
a 5-channel receiver, operating in Precision Positioning Service (PPS) mode at a rate of 1 Hz. The GPS 
data required little preprocessing, aside from extracting the data from the telemetry stream. There were 
several issues with the data itself, however. First, the east velocity channel saturated at approximately 28 
sec. after the HXLV was dropped from the B-52 carrier vehicle, due to inadequate telemetry system bit 
assignments. The east velocity data unsaturated during the vehicle descent, at a time of approximately 448 
sec. after drop. Secondly, the GPS antenna carried on the HXRV itself was activated at the instant of the 
drop (previous GPS points from the HXRV telemetry stream were provided by the vehicle flight management 
computer, but were derived from a GPS antenna carried on the B-52), and as a result, suspected blockage 
or multi path effects seem to have corrupted the first several data points recorded after the drop. Therefore, 
the first several points were deleted from the data set and were not considered in the reconstruction process. 
Also, there were several brief dropouts of all GPS data streams along the trajectory, as well as several stale 
frames and repeated time points during the descent, all of which were deleted from the data set. 

The GPS measurement noise characteristics were taken from Ref. 54 and were scaled linearly by the 
Figure of Merit (FOM), which is a telemetry parameter that takes on a value depending on the estimated 
errors in the GPS-based navigation solution. The FOM depends on many factors including the number of 
satellites used to compute the navigation solution and the dilution of precision in the event of an unfavorable 
geometry between the satellites and the receiver. 

3. Radar Data 

Many radar stations were tasked to provide measurements of the Hyper-X vehicle during its flight both for 
range safety and for trajectory reconstruction purposes. This tracking data is summarized in Ref. 55 along 
with details of the corrections made for atmospheric refraction. The tracking data must be further corrected 
for unknown bias and scale factor errors in order to be useful for trajectory reconstruction. 

During the captive carry portion of the flight, highly accurate GPS-aided INS navigation solutions are 
available in terms of geodetic altitude, geodetic latitude, and longitude. The radar tracking data are also 
available during the same portion of the flight, providing position measurements of the target relative to 
the tracking station, in terms of range, azimuth, and elevation measurements. The GPS/INS data may be 
related to the radar tracking data by means of the transformations given in Sec. IV D. A simple systematic 
error model relating measured quantities (v m , n m , and Xm) to the true quantities (v, k, and \) for the 
tracking data is of the form 


(1 + tv) V + Vv 

(55) 

(1 + tn) « + 

(56) 

(1 + t X ) X + V X 

(57) 


where £ K , and are the scale factor errors associated with the range, azimuth, and elevation data, and 
rj u , r) K , and r/ x are bias errors. By defining the radar error parameter state space as 


x r 


{ tv t. tx 


Vv Vk Vx 


(58) 


the model for the measurement equation in terms of the navigation data may be written as 
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Figure 7. Radar and GPS Residuals as Seen by MOTR-16 


where v nav , K nav , and Xnav are the range, azimuth, and elevation as computed from the navigation solution, 
w r = | w Vnav w Knav w Xriav w v w K w x | represents the zero-mean uncertainty associated with the 
range, azimuth, and elevation as computed from the INS data (essentially the navigation uncertainty trans- 
formed into range, azimuth, and elevation coordinates) and the measurement noise uncertainty for the radar 
measurements themselves. These uncertainties may be characterized in terms of their covariance matrix, 
R, r = E (w r w^) . 

The error parameters were estimated by using a Kalman filter implementation, as is described in Sec. II A. 
(A similar method for radar calibration using GPS measurements was used in Ref. 56 for air traffic control 
applications.) The resulting minimum-variance estimates of the error parameters were then applied to the 
rest of the radar tracking data in order to provide calibrated measurements of the vehicle trajectory during 
the flight experiment. Unfortunately, the tracking data provided were at low elevation angles and as a result, 
suffered from high noise content, atmospheric, and nonlinear effects that could not be accounted for in the 
calibration method. Therefore, the radar data was deemed to be unsuitable for direct use in the trajectory 
reconstruction. The data did prove useful, however, for providing reasonableness checks for the GPS data. 

As an example of the reasonableness checks, consider the MOTR-16 tracking station, located at Van- 
denberg Air Force Base (VAFB). This station provided tracking measurements on a C-Band transponder 
located on the starboard side of the HXRV. The GPS measurements of altitude, latitude, and longitude may 
be transformed into range, azimuth, and elevation data, as would be seen from this tracking station. Fig. 7 
shows the residuals between the transformed GPS data and the calibrated radar measurements. Although 
the range and elevation residuals seem to have some systematic errors (i.e., nonlinearity), and results of the 
comparison are noisy, the plots serve to show that the residuals between the GPS-derivecl quantities and the 
radar measurements do not seem to change in character between periods of relatively benign flight, such as 
that seen by the Hyper-X vehicle shortly after drop, and periods of high Mach number and high acceleration 
such as that encountered toward the end of the boost segment. Therefore, these results serve to give some 
credibility to the use of the GPS data for trajectory reconstruction purposes. 
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Longitude, deg 


Figure 8. Synoptic Weather Model 


B. Atmosphere Reconstruction 

Unlike the planet relative trajectory reconstruction, a typical atmosphere reconstruction relies solely upon 
measurements taken with large space and time separations between the measurement point and the vehicle 
during its flight. This fact is no different for the case of the Hyper-X flight. A total of 16 weather balloons 
were launched in order to provide the atmospheric data required to compute parameters of interest that are 
not estimated in the trajectory reconstruction (since the methods reconstruct a planet-relative trajectory), 
such as the dynamic pressure, Mach number, and so on. A total of 4 balloons were launched from each of 
Edwards Air Force Base (EAFB), VAFB, Point Magu Naval Air Weapons Center, and San Nicolas Island 
Naval Station. The map in Fig. 8 shows the separation of the balloon launch sites from the HXLV drop 
site. The balloon data was assembled into a synoptic model containing the ambient atmospheric properties 
as a function of altitude above the point <j> g d TBj = 33.45°, 9 re f = —120.0°, which is marked by the X on 
the map in Fig. 8. The HXLV stack drop and HXLV/HXRV separation points are also shown for reference. 
Weather data at a point other than the location of the synoptic model reference point must be transformed 
by making use of pressure gradient information to map the pressure to such a location. The temperature 
table is assumed to be invariant with respect to latitude and longitude, and the density follows directly from 
the application of the equation of state. 

C. Results of the Trajectory Reconstruction 

The results of the application of the EKF method to the Hyper-X trajectory reconstruction are shown in the 
following Figs. 9- 12. Fig. 9 shows the north velocity during the boost phase (ending at the instant before 
the opening of the engine cowl door) as a function of time since drop, which occurred at 50396.686 sec since 
midnight Pacific Standard Time (PST) on March 27, 2004. The plot shows the EKF-derived solution with 
a solid curve, the on-board navigation solution (with some prost-processing, as discussed in Ref. 58) with a 
dashed curve, and the GPS velocity measurement. In this plot, the IMU drift is clear, and the improvement 
made by using the EKF method can be seen after approximately 50 sec after drop. 

The improvement can be seen more clearly in Fig. 10, which shows the residuals between the GPS 
measurement and the EKF solution in the solid curve with circles at the data points, and the residuals 
between the GPS measurement and the on-board solution with the solid curve with squares at the data 
points. The dashed curves show the +/- one standard deviation levels in the north velocity as output from 
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Figure 9. North Velocity Comparison During Boost 


the EKF method. 

Similar results occur for the other components of velocity and position measurements. The improvements 
are summarized in Table 3, which shows the Root Mean Square (RMS) errors between the EKF solution 
and the GPS solution, and the same with the on-board navigation solution. An improvement is made by 
the application of the EKF method in all measurement components, although some cases are more dramatic 
than others. The EKF-derived component to show the smallest improvement over the on-board navigation 
solution is the east velocity. This result is perhaps not surprising, since the east velocity measurement 
channel saturated during the boost, which limited the number of useful data points to be processed. 


Table 3. Summary of RMS Errors 


Measurement 

EKF 

NAV 

Ratio 

Observations 

h g d 

21.59 ft 

42.69 ft 

1.977 

571 

4>gd 

1.202-lCT 5 ° 

3.788-10" 4 ° 

31.51 

571 

e 

2.346-10" 4 ° 

5.779-10 -4 ° 

2.463 

571 

u gd 

0.08972 ft/s 

0.3517 ft/s 

3.920 

571 

v gd 

0.9515 ft/s 

0.9746 ft/s 

1.024 

164 

W gd 

0.3189 ft/s 

0.5893 ft/s 

1.848 

571 


Finally, the time histories of the atmosphere-relative Mach number, angle of attack, and dynamic pressure 
are shown in Fig. 11 for the boost portion of the flight and Fig. 12 for the descent phase of the flight, which 
began at the instant after the engine cowl door closed. Data from the engine test phase itself is sensitive 
and cannot be shown. The various parameter identification maneuvers during the descent can be clearly 
seen in the angle of attack time history. The plots show the estimates for these quantities with a solid curve, 
and the +/- 3 standard deviation values in the dashed curves. It should be noted that these uncertainties 
are coupled in terms of the altitude uncertainty, and therefore the estimated Mach number uncertainty, 
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Figure 10. North Velocity Residual During Boost 


for example, includes the planet-relative velocity uncertainty, wind uncertainty for a specific altitude, wind 
uncertainty due to altitude uncertainty, speed of sound uncertainty for the specific altitude, and speed of 
sound uncertainty due to uncertainty in the altitude. All of these factors were taken into consideration when 
computing the standard deviations of the trajectory. 

VI. Conclusions 

This paper has discussed the formulation and implementation of an Extended Kalman Filtering algorithm 
for use in reconstructing the flight path of the Hyper-X/X-43A vehicle, using in-flight inertial measurements 
in the dynamic model combined with redundant measurements of the vehicle state from independent sources, 
particularly from the Global Positioning System. The equations of motion were formulated in a local, non- 
inertial coordinate system, and included the effects of several systematic error sources. Equations relating 
the state to measured quantities were provided for a variety of data sources of interest to this problem. Data 
recorded during the Hyper-X flight 2 on March 27, 2004, were first preprocessed and then processed using the 
estimation methods discussed in this paper, leading to an improvement over the on-board navigation solution. 
Results for the boost and descent phases of the flight were shown. The trajectory estimate determined using 
the EKF method was chosen as the project Best Estimated Trajectory (Ref. 59), principally because the 
method allowed for a systematic procedure with which to consider all sources of uncertainty in the various 
measurements, and as a result was able to provide the uncertainty values for the trajectory itself. 
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Figure 11. Mach Number, Angle of Attack, and Dynamic Pressure During Boost 
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